\section{基于深度学习特征提取匹配的全景三维重建算法设计与实现}\label{algorithm}


本章将详细介绍本文三维重建系统的实现，包括基于Superpoint的特征点提取、基于Superglue的特征点匹配、基于Colmap的SFM与MVS。
\subsection{概述}
由于Colmap在研究者与开发人员中保持着相当的热度，Colmap格式的数据可以接入到丰富的算法框架中，为测试与算法设计提供了很大的便利，
因此本文在Colmap的基础上，通过对其接口的复用，实现了一个多模块可替换的三维重建工具箱，如图\ref{fig:pipeline}。在多种算法的比较下
我们发现结合深度学习的模块对室内场景的重建有一定的提升，因此接下来将主要介绍最优的流程。

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/pipeline.png}
    \caption{三维重建流程}
    \label{fig:pipeline}
\end{figure}

我们在具体实现的过程中利用了面相对象的编程思想。对于特征提取匹配模块我们继承了pytorch的nn.Model类，构成BaseModel抽象类。所有的特征提取、匹配算法需要继承自该类并实现两个接口：forward与init。这样的好处在于
我们可以使用一套代码动态地进行多种算法的组合。所有算法模块的输出转换为Colmap的格式进行SFM。由于Colmap受到大量MVS算法支持，使得我们可以接入到各种MVS框架中，
从而实现整体各个模块可替换的效果。BaseModel代码如下：
\newpage
\begin{lstlisting}
    class BaseModel(nn.Module, metaclass=ABCMeta):
        default_conf = {}
        required_inputs = []
        def __init__(self, conf):
            """调用子类的初始化函数"""
            super().__init__()
            self.conf = conf = {**self.default_conf, **conf}
            self.required_inputs = copy(self.required_inputs)
            self._init(conf)
        def forward(self, data):
            """检查数据并调用子类的forward函数"""
            for key in self.required_inputs:
                assert key in data, 'Missing key {} in data'.format(key)
            return self._forward(data)
        @abstractmethod
        def _init(self, conf):
            """子类需要实现的抽象接口"""
            raise NotImplementedError
        @abstractmethod
        def _forward(self, data):
            """子类需要实现的抽象接口"""
            raise NotImplementedError
\end{lstlisting}\label{list:basmodel}
\subsection{基于Superpoint实现特征点提取}

Superpoint\cite[]{detone2018superpoint}是自监督特征点检测网络。其具有较好的鲁棒性，在弱纹理区域也具有较好的表现。由于其自监督的特性，可以进行大规模的训练而不需要提供真值。

\subsubsection{网络架构}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/superpoint.png}
    \caption{Superpoint网络框架}
    \label{fig:sp}
\end{figure}
网络架构如图\ref{fig:sp}。
Superpoint采用编码-解码结构。输入图像经过共享编码器之后输入到特征点解码器得到特征点概率图，输入到描述子解码器后得到描述子图。

\textbf{共享编码器：}共享编码器由三组卷积神经网络构成。每一层包括卷积、激活、池化，最终得到（$\frac{H}{8}$，$\frac{W}{8}$，C）大小的特征图。

\textbf{特征点解码器：} 将共享编码器输出的特征图再进行两次卷积、激活，最终变成（65,$\frac{H}{8}$，$\frac{W}{8}$）的张量。其中65个channel的最后一个为dustbin，
用来去除较差的结果。

\textbf{描述子解码器：} 将共享编码器输出的特征图再进行两次卷积、激活，并且对其进行归一化得到（256,$\frac{H}{8}$，$\frac{W}{8}$）大小的归一化张量图，作为描述子图。

\textbf{损失函数：} 损失函数包括特征点loss与描述子loss。
\begin{itemize}
    \item \textbf{特征点损失函数：}首先将Label进行降采样，分解成$\frac{H}{8}$、$\frac{W}{8}$的大小，接着对网络输出的特征点图进行Softmax，最后按照公式\ref*{equ:cel}计算平均交叉熵。
    \begin{equation}
        \label{equ:cel}
        \begin{split}
            l_{p}(x_{hw},y)&=-\log(\frac{\exp(x_{hwy})}{\sum^{65}_{k=1}\exp(x_{hwk})})\\
            L_{p}(X,Y)&=\frac{1}{H_{c}W_{c}} \sum^{H_{c}W_{c}}_{h=1,w=1} l_{p}(x_{hw};y_{hw})
        \end{split}
    \end{equation}
    
    这里$H_{c}W_{c}$为网络输出的特征图。y为降采样之后的Label。
    \item \textbf{描述子损失函数：} 我们希望描述子可以在图像经过变换之后对于同一特征点仍然有相似的表现。因此描述子损失函数是对源图像描述子图与变换后描述子图作用的函数。
    首先定义了一个对应集用来判断源图像与变换图像像素点的对应关系：
    \begin{equation}\label{equ:kptloss}
        s_{hwh'w'}=\left\{
            \begin{array}{lr}
               1,\Vert H_{P_{hw}}-p_{h'w'} \Vert \leq 8\\
               0,\Vert H_{P_{hw}}-p_{h'w'} \Vert > 8
               \end{array}
            \right.
    \end{equation}
    上式含义为两个像素点的欧氏距离小于等于8认为这两个像素是对应的。接下来对于源图像经过网络输出的描述子图的像素与变换图经过网络输出的
    描述子图的每一个像素有：
    \begin{equation}\label{equ:descloss}
        \begin{split}
            l_{d}(d,d';s)&=\lambda_{d}*s*max(0,m_{p}-d^{T}d')+(1-s)*max(0,d^{T}d'-m_{n})\\
            L_{d}(D,D',S)&=\frac{1}{(H_{c}W_{c})^{2}}\sum_{h=1,w=1}^{H_{c}W_{c}} \sum_{h'=1,w'=1}^{H_{c}W_{c}} l_{d}(d_{hw},d'_{h'w'};s_{hwh'w'})
        \end{split}
    \end{equation}
    描述子Loss可以分为正对应loss与负对应loss，当两个描述子对应时希望他们的内积尽可能小，而当其不对应时希望他们的内积尽可能大。公式中有三个参数，其中$\lambda_{d}$用于调节正对应与负对应损失的平衡关系，
    $m_{n}$为负对应边界，$m_{p}$为正对应边界。
\end{itemize}

最后综合公式\ref{equ:descloss}与公式\ref{equ:kptloss}，整体损失函数为：
\begin{equation}
    L(X,X',D,D';Y,Y',S)=L_{p}(X,Y)+L_{p}(X',Y')+\lambda L_{d}(D,D',S)
\end{equation}




\subsubsection{训练流程}


\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/训练流程.png}
    \caption{Superpoint训练流程}
    \label{fig:train_pipeline}
\end{figure}

\textbf{人工数据集预训练：} Superpoint首先通过程序构建各种简单的多边形并为其配置标签，输入到网络中进行训练，从而得到一个初步预测特征点的
网络。

\textbf{真实数据集增强：} 在输入真实数据集之前，每张图片都随机生成了100个单应性矩阵，从而产生了100张变换之后的图片，这大大增强了我们的数据集。同时使得
检测到的特征点具有更好的鲁棒性。

\textbf{无监督学习：} 将源图像与变换图一起输入到网络中，将预先在人工训练集上训练好的网络预测的结果作为伪真值，进行训练。

\subsubsection{代码实现}

在实现过程中我们用pytorch重新实现了Superpoint的模型并进行了训练，完成了复现。在此基础上首先我们构造一个继承自BaseModel的类，接着将init函数定义为Superpoint的初始化，主要包括关键点判决阈值、最大关键点数量、边界大小。将forward函数定义为网络进行一次前向传播的过程。返回值包括Keypoints(特征点坐标)、Describers（每个特征点的描述子）、scores（特征点质量）。
在实际进行预测的时候，只需要动态实例化一个Model类，将superpoint的配置导入进其中，输入图片，即可得到该图片的特征点。SIFT等其他特征点也可以按照同样的方式实现。

\textbf{子类：}

\begin{lstlisting}
    class SuperPoint(BaseModel):
        default_conf = {
            'nms_radius': 4,
            'keypoint_threshold': 0.005,
            'max_keypoints': -1,
            'remove_borders': 4,
        }
        required_inputs = ['image']
        def _init(self, conf):
            self.net = SP(conf)
        def _forward(self, data):
            return self.net(data)
\end{lstlisting}

\textbf{模型关键代码：}
\begin{lstlisting}
    # 共享编码器
    x = self.relu(self.conv1a(x))
    x = self.relu(self.conv1b(x))
    x = self.pool(x)
    x = self.relu(self.conv2a(x))
    x = self.relu(self.conv2b(x))
    x = self.pool(x)
    x = self.relu(self.conv3a(x))
    x = self.relu(self.conv3b(x))
    x = self.pool(x)
    x = self.relu(self.conv4a(x))
    x = self.relu(self.conv4b(x))
    # 检测子
    cPa = self.relu(self.convPa(x))
    semi = self.convPb(cPa)
    #
    prob = self.softmax(semi)
    prob = prob[:, :-1, :, :]  # remove dustbin,[B,64,H,W]
    # Reshape to get full resolution heatmap.
    prob = pixel_shuffle(prob, self.grid_size)  # [B,1,H*8,W*8]
    prob = prob.squeeze(dim=1)#[B,H,W]

    # 描述子
    cDa = self.relu(self.convDa(x))
    out = self.convDb(cDa)
    dn = torch.norm(out, p=2, dim=1)  # Compute the norm.
    desc_raw = out.div(torch.unsqueeze(dn, 1))  # Divide by norm to normalize.
    ##
    # # 描述子插值
    desc = F.interpolate(desc_raw, scale_factor=self.grid_size, mode='bilinear', align_corners=False)
    desc = F.normalize(desc, p=2, dim=1)  # normalize by channel
    prob = {'logits':semi, 'prob':prob}
    desc = {'desc_raw':desc_raw, 'desc':desc}
    return prob, desc
\end{lstlisting}

通过动态实例化模型，将数据输入即可获得特征点与描述子。

\subsection{基于Superglue实现特征点匹配}
SuperGlue\cite[]{sarlin2020superglue}是Magicleap团队在Superpoint\cite[]{detone2018superpoint}之后的续作。其目标是构建一个分配矩阵P使得所有的对应关系都有一个置信度，即输入图像对分别由M，N个特征，则对$P\in [0,1]^{M\times N}$有:

\begin{equation}
    P 1_{N}\leq 1_{M} ,P^{T} 1_{M}\leq 1_{N}
\end{equation}

Superglue对于匹配问题，将其抽象为从一个概率分布转换为另一个概率分布。分配矩阵P表示了这种转换的代价，即匹配概率。为了得到这个匹配矩阵，可以从人类视觉匹配的角度出发，人类在寻找匹配时会考虑到点的位置与像素信息、点与另一个点的相对位置关系。
另外，在另一张图中的信息也可以用于判别一些难以区分的情况。同时，我们往往会在两张图片中来回对比，将注意力集中在某些点与点的关系上。这同样也是Transformer的思想。因此根据以上的思想，Superglue采用了一种基于注意力的特征匹配网络。具体网络架构如图\ref{fig:sg_pipeline}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/SuperGluePipeline.png}
    \caption{SuperGlue流程}
    \label{fig:sg_pipeline}
\end{figure}
\subsubsection{注意力机制图神经网络}
\textbf{特征点编码：}相比于传统匹配算法， SuperGlue同时参考了特征点的位置信息与像素信息，将输入的描述子（像素信息）与特征点（位置信息）通过MLP进行编码（如公式\ref{equ:encoder}），更加合理。

\begin{equation}\label{equ:encoder}
    x_{i}=d_{i}+MLP_{enc}(p_{i})
\end{equation}

\textbf{多重图神经网络：}如图\ref{fig:sg_pipeline}所示，图中的节点为两个图像的特征点，边表示特征点之间的连接关系。存在有两类连接关系：$\varepsilon_{cross}$ 表示图像之间的连接，$\varepsilon_{self}$表示图像内的连接。
多重图神经网络从低层网络开始，对每一个节点通过同时聚合所有节点的所有给定边上的消息，在每一层计算更新的节点的状态，从而完成节点之间信息的传递。
我们用$^{l}x_{i}^{A}$表示A图下第l层的节点i，$m_{\varepsilon  \rightarrow  i}$表示所有连接i的节点聚合的信息，包括自相关与互相关。那么每一层的信息更新为：

\begin{equation}
    ^{(l+1)}x_{i}^{A}=^{(l)}x_{i}^{A}+MLP([^{l}x_{i}^{A}||m_{\varepsilon \rightarrow i}])
\end{equation}

这里$[ \cdot || \cdot]$表示连接。同时另一张图也进行了同样的更新。多重图神经网络在奇数层更新自注意力，在偶数层更新互注意力。这一过程即是接下来要介绍的注意力聚合。

\textbf{注意力聚合：}我们可以定义一种查询关系：通过键q根据属性k检索值v，则代价聚合公式：

\begin{equation}\label{equ:agg}
    m_{\varepsilon \rightarrow i } = \sum_{j:(i,j)\in \varepsilon} \alpha_{ij}v_{j}
\end{equation}

公式\ref{equ:agg}中，$\alpha_{ij}$是对于所有的键-值-属性的相似度进行softmax运算，具体计算方式如下：

\begin{equation}
    \alpha_{ij}=Softmax_{j}(q_{i}^{T}k_{j})
\end{equation}

k,q,v的计算方式是通过对特征x进行线性映射得到。假设查询点i在图片Q，源点在图像S，则有：
\begin{equation}
    \begin{split}
        q_{i}&=W_{1}^{l}x_{i}^{Q}+b_{1}\\
        \begin{pmatrix}
            k_{j}\\
            v_{j}
        \end{pmatrix}&=\begin{pmatrix}
            W_{2}\\
            W_{3}
        \end{pmatrix}^{(l)}x_{j}^{S}+\begin{pmatrix}
            b_{2}\\
            b_{3}
        \end{pmatrix}
    \end{split}
\end{equation}

每一层网络都有各自的[k,q,v]参数，对输入该层的所有特征都共享同一组参数。

\textbf{匹配代价：}
最后对每一个特征点进行线性映射得到匹配代价：

\begin{equation}
    f_{i}^{A}=W\cdot^{(L)}x_{i}^{A}+b,\forall i \in A
\end{equation}

\subsubsection{优化匹配层}
这一层主要解决从特征点到$P_{M\times N}$分配矩阵的映射问题。

\textbf{预测分数:} 预测分数通过对两个特征向量进行内积运算获得。

\begin{equation}
    S_{i,j}=<f_{i}^{A},f_{j}^{B}>,\forall (i,j)\in A \times B
\end{equation}

这里的S与f没有进行归一化以反映匹配的置信度。

\textbf{遮挡与可见性：}由于图像中并非所有特征都可以找到匹配的对象，当存在遮挡等因素会出现特征无匹配的情况，为了解决这个问题，
在P外层加了一层dustbin，所有没有匹配的点都将划分到该区域。

\textbf{并行优化：}最后，通过Sinkhorn算法进行并行计算以使模型应用到GPU上。

由于Superglue本身使用pytorch实现，因此这里直接使用了作者的模型与权重进行调用。在本文中调用方式与调用Superpoint类似。

\subsection{基于Colmap SFM实现稀疏重建}
\subsubsection{Colmap SFM原理与流程}
% 其流程如图\ref{fig:colmap_sfm_pipeline}
Colmap的SFM采用的是\ref{sec:sfm}节提到的增量式SFM。Colmap读取输入的特征点与匹配对以生成场景图（Scene Graph）和匹配矩阵（Matching 
Matrix）。之后采用增量式的 SfM 算法逐步增加视角迭代地优化重投影误差，确定不同视图和点
云间的可视关系，并得到场景中的相机位姿和表示场景结构的稀疏点云。

增量式SfM 采用迭代的算法逐步优化稀疏重建结果，其算法流程如下：
\begin{enumerate}[leftmargin=4em]
    \item \textbf{对无序图像进行特征匹配和三角测量以进行初始化。}由于场景中的点不能被持续观测，所以特征追踪可能因此导致错误匹配，进而导致后续的三
    角化发生错误，因此 Colmap 使用 RANSAC 对多帧观测的图像进行三角化。在Colmap中选择使用DLT三角化法\cite[]{hartley2003multiple}。最终点要求角度足够且深度非负。
    \item \textbf{通过已有点云重新估计相机位姿}
    \item \textbf{进行全局和局部的BA优化}BA优化的思路是通过将特征点进行重投影，以平方差的方式表示。对平方差求解偏导数或者说雅可比矩阵。
    将变量沿着雅可比矩阵方向下降，从而减小整体误差。优化对象可以是相机参数、特征点坐标等。在进行雅可比矩阵求导的时候，对于旋转矩阵的求导需要使用李代数来进行。在优化过后相机参数发生变化，因此需要重新三角化。
    \item \textbf{增量添加视角，再次进行三角测量和相机位姿估计。}
    \item \textbf{再进行BA优化修正结果}
\end{enumerate}

% \begin{figure}[H]
%     \centering
%     \includegraphics[width=0.8\textwidth]{figures/colmap_pipeline.png}
%     \caption{Colmap流程}
%     \label{fig:colmap_sfm_pipeline}
% \end{figure}

\subsubsection{colmap的接口设计}
接下来的部分主要解决特征点和匹配对如何接入Colmap 的SFM Pipeline，以及Colmap 的SFM结果与各种MVS算法的接入问题。colmap的数据管理是基于数据库进行的。所有图片、特征点、匹配对都会在数据库中进行注册、操作和管理，这种方式
具有很好的拓展性与操作性。在代码编写过程中我们参考了colmap仓库提供的python脚本，将特征点、匹配对输入到数据库中，根据其格式要求编写matchlist完成几何矫正。
colmap 提供的接口为：

\begin{itemize}[leftmargin=4em]
    \item \textbf{feature\_importer:} 特征点导入
    \item \textbf{matches\_importer:} 匹配对导入
    \item \textbf{mapper:} SFM接口。
\end{itemize}


在实际使用中我们只使用到了mapper这一接口。因为colmap提供的feature importer、matches importer对格式有一定要求，而其底层原理可以通过python与
数据库交互实现，因此我们选择用python实现这一过程。通过python调用mapper接口，将整个流程统一到python环境下便于管理。具体代码如下：

\textbf{SFM关键代码}

\begin{lstlisting}
    create_empty_db(database)
    import_images(
        colmap_path, sfm_dir, image_dir, database, single_camera)
    image_ids = get_image_ids(database)
    import_features(image_ids, database, features)
    import_matches(image_ids, database, pairs, matches,
                   min_match_score, skip_geometric_verification)
    if not skip_geometric_verification:
        geometric_verification(colmap_path, database, pairs)
    stats = run_reconstruction(
        colmap_path, sfm_dir, database, image_dir, min_num_matches)
\end{lstlisting}

\textbf{匹配对导入}
\begin{lstlisting}
    def import_matches(image_ids, database_path, pairs_path, matches_path,
                   min_match_score=None, skip_geometric_verification=False):
        logging.info('Importing matches into the database...')
        with open(str(pairs_path), 'r') as f:
            pairs = json.load(f)
        hfile = h5py.File(str(matches_path), 'r')
        db = COLMAPDatabase.connect(database_path)
        # print(pairs)        matched = set()
        for name0, name1 in tqdm(pairs):
            id0, id1 = image_ids[name0], image_ids[name1]
            if len({(id0, id1), (id1, id0)} & matched) > 0:
                continue
            pair = name0+' '+name1
            if pair not in hfile:
                raise ValueError(
                    f'Could not find pair {(name0, name1)}... '
                    'Maybe you matched with a different list of pairs? '
                    f'Reverse in file: {names_to_pair(name0, name1) in hfile}.')
            matches = hfile[pair]['matches0'].__array__()
            valid = matches > -1
            if min_match_score:
                scores = hfile[pair]['matching_scores0'].__array__()
                valid = valid & (scores > min_match_score)
            matches = np.stack([np.where(valid)[0], matches[valid]], -1)
            db.add_matches(id0, id1, matches)
            matched |= {(id0, id1), (id1, id0)}
            if skip_geometric_verification:
                db.add_two_view_geometry(id0, id1, matches)
        hfile.close()
        db.commit()
        db.close()
\end{lstlisting}

\subsection{基于Colmap MVS实现稠密重建}

Colmap的MVS主要进行了深度图计算与点云融合。其原理与\ref{sec:mvs}节接近，流程如图\ref{fig:colmap_mvs_pipeline}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{figures/colmap_mvs_pipeline.png}
    \caption{Colmap流程}
    \label{fig:colmap_mvs_pipeline}
\end{figure}

深度图估计主要包括以下几个步骤：

\begin{enumerate}[leftmargin=4em]
    \item  筛选最优的参考图像，找出视差计算的选集
    \item 计算参考图像上其他图像的匹配代价
    \item 将代价进行聚合计算深度值
    \item 滤波以优化深度图
\end{enumerate}

Colmap在MVS方面采用的是PatchMatchStereo（参考\ref{sec:mvs_pms}节），使用NCC（归一化互相关/Normalized Cross Correlation）作为代价计算的指标，通过PatchMatch进行信息传递。
其后处理方法是利用光学一致性与几何一致性约束对深度图进行过滤。最后对深度图进行融合得到稠密点云。相比于原始PMS，Colmap在最优帧选取上不仅考虑到了三角化点的数量，也考虑到了点的
分布情况对精度的影响。这一点在后续实验中也会得到验证。

在实际操作中，我们选择使用python调用Colmap的接口，包括以下：

\begin{itemize}[leftmargin=4em]
    \item \textbf{image\_undistorter：} 图像去畸变
    \item \textbf{patch\_match\_stereo：} PMS计算深度图
    \item \textbf{stereo\_fusion：} 深度图融合
\end{itemize}

通过python调用上述接口可以实现Colmap的MVS，也可以通过调用openmvs的命令行实现openmvs的MVS。


\subsection{本章小结}

本章介绍了本文选取的方法。首先介绍了本文选取的深度学习特征点提取与匹配算法的网络架构，
接着介绍了深度学习匹配对、特征点到Colmap数据格式的转换。接着介绍了如何使用Colmap实现SFM与MVS。关键部分代码在文中有展示，具体代码请参考代码仓库。